loo

libraries

library(tidyverse)
library(brms)
library(bayesplot)
library(here)
library(glue)
library(scales)

theme_set(theme_light())

source(glue("{params$common_dir_str}/brms_model.R"))

data

obs_only <- 
  read_csv(glue("{params$model_dir_str}/data/stimulation_obvs.csv")) %>%
  mutate(subj = as_factor(subj),
         obs_degree = error,
         error = obs_degree * (pi/180))

## Parsed with column specification:
## cols(
##   subj = col_double(),
##   subj_index = col_double(),
##   stimulation = col_double(),
##   error = col_double()
## )

full model

# define prior

bprior <- 
  prior(normal(4, 0.5), class = "b", coef = "intercept", dpar = "circSD") + 
  prior(normal(0, 0.5), class = "b", coef = "stimulation", dpar = "circSD") +
  prior(normal(0, 0.5), class = "sd", coef = "Intercept", dpar = "circSD", group = "subj_index") + 
  prior(normal(0, 0.5), class = "sd", coef = "stimulation", dpar = "circSD", group = "subj_index") + 
  prior(normal(0, 1.5), class = "b", coef = "intercept", dpar = "theta") + 
  prior(normal(0, 1), class = "b", coef = "stimulation", dpar = "theta") +
  prior(normal(0, 1.5), class = "sd", coef = "Intercept", dpar = "theta", group = "subj_index") + 
  prior(normal(0, 1), class = "sd", coef = "stimulation", dpar = "theta", group = "subj_index")


iter = 4000
warmup = 2000
cores = 4
chains = 4
n_post_samples = (iter - warmup) * chains

fit_full <- brm(bform_full, obs_only, family = vm_uniform_mix, prior = bprior, stanvars = stanvars,
                 warmup = warmup, iter = iter, cores = cores, chains = chains, 
                 control = list(adapt_delta = 0.99), inits = 0, 
                 file =  glue("{params$save_dir_str}/fit_full"))

fit check

divergences

model_fit <- fit_full

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))

##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

reduced model: no pre to post change on circSD parameter

# define prior

bprior <- 
  prior(normal(4, 0.5), class = "b", coef = "intercept", dpar = "circSD") + 
  #prior(normal(0, 0.5), class = "b", coef = "stimulation", dpar = "circSD") +
  prior(normal(0, 0.5), class = "sd", coef = "Intercept", dpar = "circSD", group = "subj_index") + 
  prior(normal(0, 0.5), class = "sd", coef = "stimulation", dpar = "circSD", group = "subj_index") + 
  prior(normal(0, 1.5), class = "b", coef = "intercept", dpar = "theta") + 
  prior(normal(0, 1), class = "b", coef = "stimulation", dpar = "theta") +
  prior(normal(0, 1.5), class = "sd", coef = "Intercept", dpar = "theta", group = "subj_index") + 
  prior(normal(0, 1), class = "sd", coef = "stimulation", dpar = "theta", group = "subj_index")


iter = 4000
warmup = 2000
cores = 4
chains = 4
n_post_samples = (iter - warmup) * chains

fit_DcircSD_null <- brm(bform_DcircSD_null, obs_only, family = vm_uniform_mix, prior = bprior, stanvars = stanvars,
                 warmup = warmup, iter = iter, cores = cores, chains = chains, 
                 control = list(adapt_delta = 0.99), inits = 0, 
                 file =  glue("{params$save_dir_str}/fit_DcircSD_null"))

fit check

divergences

model_fit <- fit_DcircSD_null

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))

##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

reduced model: no pre to post change on pMem parameter

# define prior

bprior <- 
  prior(normal(4, 0.5), class = "b", coef = "intercept", dpar = "circSD") + 
  prior(normal(0, 0.5), class = "b", coef = "stimulation", dpar = "circSD") +
  prior(normal(0, 0.5), class = "sd", coef = "Intercept", dpar = "circSD", group = "subj_index") + 
  prior(normal(0, 0.5), class = "sd", coef = "stimulation", dpar = "circSD", group = "subj_index") + 
  prior(normal(0, 1.5), class = "b", coef = "intercept", dpar = "theta") + 
  #prior(normal(0, 1), class = "b", coef = "stimulation", dpar = "theta") +
  prior(normal(0, 1.5), class = "sd", coef = "Intercept", dpar = "theta", group = "subj_index") + 
  prior(normal(0, 1), class = "sd", coef = "stimulation", dpar = "theta", group = "subj_index")


iter = 4000
warmup = 2000
cores = 4
chains = 4
n_post_samples = (iter - warmup) * chains

fit_DpMem_null <- brm(bform_DpMem_null, obs_only, family = vm_uniform_mix, prior = bprior, stanvars = stanvars,
                 warmup = warmup, iter = iter, cores = cores, chains = chains, 
                 control = list(adapt_delta = 0.99), inits = 0, 
                 file = glue("{params$save_dir_str}/fit_DpMem_null"))

fit check

divergences

model_fit <- fit_DpMem_null

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))

##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

reduced model: no pre to post change on pMem or circSD parameter

# define prior

bprior <- 
  prior(normal(4, 0.5), class = "b", coef = "intercept", dpar = "circSD") + 
  #prior(normal(0, 0.5), class = "b", coef = "stimulation", dpar = "circSD") +
  prior(normal(0, 0.5), class = "sd", coef = "Intercept", dpar = "circSD", group = "subj_index") + 
  prior(normal(0, 0.5), class = "sd", coef = "stimulation", dpar = "circSD", group = "subj_index") + 
  prior(normal(0, 1.5), class = "b", coef = "intercept", dpar = "theta") + 
  #prior(normal(0, 1), class = "b", coef = "stimulation", dpar = "theta") +
  prior(normal(0, 1.5), class = "sd", coef = "Intercept", dpar = "theta", group = "subj_index") + 
  prior(normal(0, 1), class = "sd", coef = "stimulation", dpar = "theta", group = "subj_index")


iter = 4000
warmup = 2000
cores = 4
chains = 4
n_post_samples = (iter - warmup) * chains

fit_DcircSD_DpMem_null <- 
              brm(bform_DcircSD_DpMem_null, obs_only, family = vm_uniform_mix, prior = bprior, stanvars = stanvars,
                 warmup = warmup, iter = iter, cores = cores, chains = chains, 
                 control = list(adapt_delta = 0.99), inits = 0, 
                 file = glue("{params$save_dir_str}/fit_DcircSD_DpMem_null"))

fit check

divergences

model_fit <- fit_DcircSD_DpMem_null

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))

##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

loo compare

#need to do this only once?
expose_functions(fit_full, vectorize = TRUE)

fit_full <- add_loo(fit_full, file = glue("{params$save_dir_str}/fit_full"))
fit_DcircSD_null <- add_loo(fit_DcircSD_null, file = glue("{params$save_dir_str}/fit_DcircSD_null"))
fit_DpMem_null <- add_loo(fit_DpMem_null, file = glue("{params$save_dir_str}/fit_DpMem_null"))
fit_DcircSD_DpMem_null <- add_loo(fit_DcircSD_DpMem_null, file = glue("{params$save_dir_str}/fit_DcircSD_DpMem_null"))

print("fit_full")

## [1] "fit_full"

fit_full$loo

## 
## Computed from 8000 by 504 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -804.7 15.9
## p_loo         7.2  0.2
## looic      1609.5 31.8
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

print("fit_DcircSD_null")

## [1] "fit_DcircSD_null"

fit_DcircSD_null$loo

## 
## Computed from 8000 by 504 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -805.2 15.8
## p_loo         7.1  0.2
## looic      1610.4 31.6
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

print("fit_DpMem_null")

## [1] "fit_DpMem_null"

fit_DpMem_null$loo

## 
## Computed from 8000 by 504 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -804.3 15.9
## p_loo         6.7  0.2
## looic      1608.5 31.8
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

print("fit_DcircSD_DpMem_null")

## [1] "fit_DcircSD_DpMem_null"

fit_DcircSD_DpMem_null$loo

## 
## Computed from 8000 by 504 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -804.7 15.8
## p_loo         6.8  0.2
## looic      1609.5 31.7
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

loo_compare(fit_full, fit_DcircSD_null, fit_DpMem_null, fit_DcircSD_DpMem_null, criterion = "loo")

##                        elpd_diff se_diff
## fit_DpMem_null          0.0       0.0   
## fit_full               -0.5       0.2   
## fit_DcircSD_DpMem_null -0.5       0.6   
## fit_DcircSD_null       -0.9       0.7